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Results from previous shell buckling studies are presented that illustrate some of 
the unique and powerful capabilities in the STAGS finite element analysis code that 
have made it an indispensable tool in structures research at NASA over the past few 
decades. In particular, prototypical results from the development and validation of 
high-fidelity buckling simulations are presented for several unstiffened thin-walled 
compression-loaded graphite-epoxy cylindrical shells along with a discussion on the 
specific methods and user-defined subroutines in STAGS that are used to carry out 
the high-fidelity simulations. These simulations accurately account for the effects of 
geometric shell-wall imperfections, shell-wall thickness variations, local shell-wall 
ply-gaps associated with the fabrication process, shell-end geometric imperfections, 
nonuniform applied end loads, and elastic boundary conditions. The analysis 
procedure uses a combination of nonlinear quasi-static and transient dynamic 
solution algorithms to predict the prebuckling and unstable collapse response 
characteristics of the cylinders. Finally, the use of high-fidelity models in the 
development of analysis-based shell-buckling knockdown (design) factors is 
demonstrated. 


I. Introduction 

T he ever-increasing need to produce lighter-weight aerospace shell structures has led to the use 
of highly optimized structural designs and advanced material systems. Buckling is often a 
primary failure mode that drives the design of light-weight aerospace shell structures such as 
cryogenic tanks and interstage structures in launch vehicles. In particular, maximum structural 
efficiency typically necessitates the use of thin highly stressed skins that can be susceptible to 
buckling failures when subjected to compression loads. Thus, buckling is an important and often 
critical consideration in the design of these structures and reliable, validated high-fidelity 
analysis tools and design criteria are needed. 

The general topic of shell buckling has garnered a tremendous amount of attention since the 
beginning of the 20 th century. From approximately 1930 to 1967, many shell-buckling 
experiments were conducted on isotropic unstiffened cylindrical shells under a variety of loading 
conditions. Typically, the experiments yielded buckling loads that were substantially lower than 
the corresponding analytical predictions, which were based on simplified linear bifurcation 
analyses of geometrically perfect shells with nominal dimensions and idealized boundary 
conditions. The pioneering works by von Karman and Tsien (1941) * 1 , by Donnell and Wan 
(1950)“ , and by Koiter (1945) identified small deviations from the idealized geometry of a 


* Senior Research Engineer, Structural Mechanics and Concepts Branch, Hampton, VA, 23681, 
USA, Senior Member AIAA. 
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shell, known as initial geometric imperfections, as a primary source of discrepancy between 
corresponding analytical predictions and experimental results. However, the analysis capabilities 
at that time could not meet the intense computational requirements of the nonlinear analyses 
needed to predict accurately the imperfection sensitivity and buckling loads of the tested shells. 
Thus, designers often used a simplified analysis procedure with empirical data to develop safe 
designs for buckling-critical structures. Specifically, the traditional approach for designing thin- 
walled buckling-resistant isotropic shell structures is to predict the buckling load of the shell with 
a linear bifurcation buckling analysis, and then to reduce this predicted load with an empirical 
“knockdown” factor. The empirical knockdown factor is intended to account for the difference 
between the predicted buckling load and the actual buckling load for the shell determined from 
tests and is often based on the “lower bound” design recommendations reported in the NASA 
SP-8007 document (1968). 4 This approach to shell buckling design remains prominent in 
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industry practice, as evidenced by the extensive use of the NASA space vehicle design criteria. ' 

From approximately 1960 to 1989, digital computers emerged and matured into a practical tool 
for research and design. With these computational tools, considerable effort was put forth to 
determine analytically the effects of initial geometric imperfections on the buckling response of 
shells. These efforts were largely based on Koiter’s asymptotic theory 3 ' 9 for shells with nominal 
dimensions and idealized boundary conditions. In addition, in-depth experimental investigations 
were conducted at institutions such as the NASA Langley Research Center and the Wright 
Aeronautical Laboratories at Wright-Patterson Air Force Base, in addition to various academic 
institutions worldwide. It was during this time period that the measurement of initial geometric 
imperfections in test specimens and high-fidelity response measurement techniques and 
apparatus were recognized as important elements in the development of design criteria for 
buckling of shells. 10 " 13 

Another significant development that occurred in the 1960s to 1980s was the development of 
computational structural analysis codes based on finite-difference or finite-element methods. 
Two prominent examples at the forefront in the development of this technology are the NASA 
Structural Analysis code (NASTRAN 14 ) and the Structural Analysis of General Shells code 
(STAGS 15 ). These codes evolved over many years and now possess robust static and dynamic 
nonlinear structural analysis capabilities, along with extensive graphics capabilities, that enable 
scientists and engineers to model and visualize experiments or actual flight hardware to an 
enormous level of detail. Prototypical examples of the application of the STAGS code can be 
seen in the in-depth, highly detailed studies of cylinder buckling response and imperfection 
sensitivity documented in Refs. 16-18 and the nonlinear and buckling analyses of the Space 
Shuttle super- light- weight external fuel tank that is documented in Refs. 19-23. 

The results from these studies showed that highly accurate predictions of the buckling response 
of thin-walled cylinders can now be achieved based on detailed nonlinear analysis and that 
analysis-based knockdown factors could become a viable replacement for the test-based design 
knockdown factors used currently. In addition, the robust nonlinear solution algorithms and 
user-defined subroutines enabled precise, systematic and efficient exploration of the 
imperfection sensitivities that were affecting the structural response. 

Based on these promising results, Hilburger and Starnes proposed an analysis-based design 
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approach in 2003 that can accurately account for the effects of initial geometric imperfections. 
With this approach, measured initial geometric imperfection data from a representative set of 
cylindrical shells are used to determine a manufacturing-process-specific imperfection signature 
for these shells. This imperfection signature is then used as input into nonlinear finite-element 
analyses. The imperfection signature represents a "first-approximation" imperfection shape that 
is suitable for developing preliminary-design data. The STAGS code, again, played a prominent 
role in developing the analysis-based design strategy. 

In recent times, this new analysis-based approach has been further developed and used to 
generate improved design knockdown factors for NASA’s next generation of heavy-lift launch 
vehicle shell structures.' These new factors have been developed for a particular class of 
integrally stiffened metallic cylinders and account for manufacturing-process-specific 
imperfection signature as well as critical design features such as longitudinal welds which have 
been shown to significantly affect the buckling response of these shells. The STAGS finite 
element code was used extensively in the path-finding work during the development of these 
new knockdown factors as well as providing benchmark results to help develop and validate 
models using other commonly used commercial finite element codes. 

The present paper is intended to highlight some of the ground-breaking analysis capabilities of 
the STAGS finite-element code that have been developed over the past few decades and that are 
enabling new design technologies such as analysis-based shell buckling knockdown factors. To 
this end, selected results from previous numerical and experimental studies on the buckling of 
unstiffened thin-walled compression-loaded graphite-epoxy cylindrical shells are presented. 
First, a set of three composite laminated cylindrical shell test articles are described. Then, a 
high-fidelity nonlinear shell-analysis procedure is describe and includes details on how various 
modeling details are included in the STAGS model including initial geometric shell-wall 
imperfections, shell-wall thickness variations, shell-end geometric imperfections, nonuniform 
applied end loads, and elastic boundary conditions. Next, selected test and analysis correlation 
data is presented to illustrate the accuracy of the high-fidelity models. Finally, the foundational 
work towards the development and implementation of new analysis-based shell buckling 
knockdown factors are described and design implications are discussed. 

II. Test Specimens, Imperfection Measurements and Test Procedure 
Test Specimens 

The test specimens discussed in this paper were fabricated from 0.005-in.-thick AS4/3502 
graphite-epoxy pre-impregnated unidirectional tape material made by Hercules, Inc. The nominal 
unidirectional lamina properties of a typical 0.005 -in. -thick ply with a fiber volume fraction of 
0.62 are as follows: longitudinal modulus Ei = 18.5 Msi, transverse modulus E 2 = 1.45 Msi, in- 
plane shear modulus G 12 = 0.813 Msi, and major Poisson’s ratio V 12 = 0.30. The material was 
laid up on a 15.75-in.-diameter mandrel and cured in an autoclave to form three shells with 
different laminates including an 8-ply axially stiff [-45/+45/02] s laminate, an 8-ply 
circumferentially stiff [-45/+45/902] s laminate, and an 8-ply quasi-isotropic [-45/+45/0/90] s 
laminate. The resulting three shells are referred to herein as shells or specimens Cl, C2, and C3, 
respectively. These specimens had a nominal length of 16.0 in. and a nominal radius of 8.0 in. 
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The specimens had a nominal shell-wall thickness of 0.04 in. and a shell-radius-to-thickness ratio 
of 200. Both ends of the specimens were potted in an aluminum-filled epoxy resin to assure that 
the ends of the specimens did not fail prematurely during the test. The potting material extended 
approximately 1 .0 inch along the length of the specimens at each end resulting in a test section 
that was approximately 14.0 in. long. The ends of the specimens were machined flat and parallel 
to assure proper load introduction during the tests. A photograph of a typical specimen is shown 
in Figure 1 and the specimen coordinate system and geometry is shown in Figure 2. The shell 
length, test-section length, radius, and thickness are designated as L, Lj. R and t, respectively. 

Imperfection Measurements 

Three-dimensional surveys of the inner and outer shell-wall surfaces of the specimens were 
made prior to testing to determine the initial geometric imperfection shapes of the shell-wall and 
shell-wall thickness distributions. Measurements were taken over a uniform grid with increments 
of 0.125 in. in the axial direction and approximately 1° of arc in the circumferential direction 
over the exposed surfaces of the specimens. The inner surface measurement was used to 
determine the initial geometric imperfection shape of a specimen, and the difference between the 
outer and inner surface measurements was used to determine the shell-wall thickness 
distribution. A contour plot of the normalized initial geometric imperfection for specimen C3 is 
shown in Fig. 3. The measured geometric imperfection w 0 is normalized by the average 
measured shell-wall thickness t ave = 0.0381 inches. These results indicate that the initial 
geometric imperfection is periodic in the circumferential direction and has slight deviations 
along the axial direction. The amplitude of the imperfection varies from +1.341f ave to -1.535W • 
A contour plot of the normalized shell-wall thickness variation for specimen C3 is shown in Fig. 
4, where the measured thickness value t 0 is normalized by the average measured shell-wall 
thickness t ave . These results indicate that the shell-wall thickness, and hence the laminate 
stiffnesses, varies significantly over a short distance. The thickness varies from 0.928 to 1.321 
times tave • Most of the thickness variation is attributed to local variations in the resin content of 
the laminate associated with the fabrication process. However, the dark blue angular pattern in 
the thickness distribution is attributed to small gaps between adjacent pieces of graphite-epoxy 
tape in some of the laminate plies that were generated during the lay-up and curing processes. 
Such a region is referred to herein as a lamina ply-gap or a ply-gap. In the particular case of shell 
C3, such a locally thin region in the shell wall consists of a 7-ply-thick laminate rather than the 
nominal 8-ply-thick laminate. For the case where one ply-gap intersects another ply-gap, the 
shell wall consists of a 6-ply-thick laminate. Moreover, these locally thin shell-wall regions have 
a significant shell-wall mid-surface eccentricity. A circumferential ply-gap is caused by a gap 
between two adjacent 90° pieces of tape, a helical ply-gap is caused by a gap between two 
adjacent 45° pieces of tape, and an axial ply-gap is caused by a gap between two adjacent 0° 
pieces of tape. Magnified cross-sectional views of typical ply-gaps in a -45° outer ply of a 
specimen laminate and in a 90° mid-surface ply are shown in Figs. 5a and 5b, respectively. The 
widths of the ply-gaps shown in Figs. 5a and 5b are on the order of 0.02 in. or approximately 
equal to half the shell-wall thickness, and the ply-gap depth is approximately 0.005 in. or 
approximately equal to the nominal lamina ply thickness. Lamina ply-gaps with gap widths as 
large as 0.1 in. have been observed in some of the shell specimens. The red colored 
longitudinally aligned patterns in the thickness contour plot are caused by locally thickened 
regions of the outermost plies of the laminate that develop during the curing process to form 
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outer shell-wall surface ridges. A magnified cross-sectional view of such a region is shown in 
Fig. 5c. Some of the shell-wall thickness features, such as lamina ply-gaps, are smaller than the 
imperfection-measurement grid spacing used in this study, and as a result, some of the smaller 
thickness variation features may not have been included in the measurements. 

Measurements of the top and bottom loading surfaces of the specimens were made every degree 
around the circumference of the specimens to determine the variation in the shell-end or loading- 
surface geometry. Typical top and bottom shell-end or loading-surface geometry variations for 
specimen C3, denoted by d top ( 6) and dbot( 6), respectively, are shown in Fig. 6. The maximum 
amplitude of this loading- surface variation is approximately 0.0015 inches, which is 
approximately 0.01% of the specimen length. 

Test Apparatus and Procedure 

The specimens were instrumented with electrical resistance strain gages, and direct-current 
differential transducers (DCDT’s) were used to measure displacements. Three non-collinear 
DCDT’s were positioned at three corners of the upper loading platen of the test machine and 
used to measure the end-shortening displacement A and the rotations (p y and (p- of the loading 
platen as illustrated in Fig. 2. Typical measured platen rotations are shown in Figure 7. To 
control the load introduction into the specimens, the upper loading platen was aligned with the 
loading surface of the specimen as well as possible before the test by adjusting leveling bolts in 
the corners of the upper platen until strains measured by selected strain gages on the specimens 
indicated a uniform axial strain distribution around the circumference of the shell. The specimens 
were loaded in compression with a 300,000-lb hydraulic universal-testing machine by applying 
an end- shortening displacement to the loading surfaces of the specimens. The specimens were 
loaded until general instability or failure of the shells occurred. 

A shadow moire interferometry technique was used to observe the shell-wall prebuckling, 
buckling and post-buckling radial deformation patterns. All data were recorded with a data 
acquisition system, and the moire patterns were recorded photographically and on videotape. 


III. Finite-Element Models and Analysis Procedure 

The shells considered in this study were analyzed with the STAGS (STructural Analysis of 
General Shells) nonlinear shell analysis code. 15 STAGS is a research-oriented finite-element 
code developed for the nonlinear static and dynamic analysis of general shells, and can include 
the effects of geometric and material nonlinearities in the analysis. It possesses a typical suite of 
elements including beams, plates, and solids, elastic and rigid links, as well as generalized 
contact elements, sandwich elements, and user-defined element routine similar to that in Abaqus. 
In addition, STAGS includes progressive damage analysis capabilities including crack 
propagation, delamination modeling via decohesion elements and property degradation methods 
for composite materials. The STAGS code also has a unique feature in that it includes several 
different user-defined subroutines that can be linked to the STAGS program executable to 
accomplish many different modeling and data processing tasks. In addition, the STAGS source 
code is open for modification by the user, which enables the implementation of new solution 
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techniques directly into the code. The user-written subroutine capability is used extensively by 
NASA researchers and in the work presented herein. 

Finite-Element Models 

A typical finite-element model of a specimen is illustrated in Fig. 2. The standard 410 
quadrilateral shell element from the STAGS element library was used in the models. The 
element is based on the Kirchoff-Love shell theory and the nonlinear Lagrangian strain tensor. 
The elements of the finite-element mesh are approximately 0.2-in. by 0.2-in. square. Each 
element has four integration points, which are distributed in such a way as to provide a modeling 
resolution of approximately 0.1-in. by 0.1-in square. This integration point spacing is on the 
order of the measurement-point spacing used to measure the initial geometric imperfections of 
the specimens. This highly refined mesh is necessary to model rapidly varying geometric and 
material parameters such as nonuniform shell-wall thicknesses and lamina stiffness properties. A 
typical finite-element model contained approximately 100,000 degrees of freedom. 

Results from geometrically perfect and imperfect shells are presented in this paper. The 
geometrically perfect shells are modeled assuming nominal shell geometry, laminate thickness, 
lamina mechanical properties, and boundary conditions. The boundary conditions consist of 
setting the circumferential and radial displacements v and w equal to zero in the l.O-in.-long 
potted boundary regions of the shell illustrated in Fig. 1, setting u(L/2, 6) = 0, and applying a 
uniform end- shortening u( -L/2, 6) = A. 

The geometrically perfect finite -element models were then modified to include the effects of the 
measured shell imperfections in order to simulate the response of the specimens more closely. 
These modeling modifications included the effects of the measured initial geometric shell-wall 
imperfections, shell-wall thickness variations, local shell-wall lamina ply-gaps, thickness- 
adjusted lamina properties, elastic boundary conditions, shell-end geometric imperfections, and 
nonuniform end loads. 

The initial geometric shell- wall imperfection w 0 (x, 6) is included in the finite-element models by 
introducing an initial radial perturbation to each node of the mesh by using the user-written 
subroutine DIMP available in STAGS. The DIMP subroutine statement is as follows (variables 
used in the subroutine are highlighted in bold font): 

subroutine DIMP(iunit, inode , prop, XYZg, XYs, vv) 

where inode is the node number from the model, XYs is a vector that contains the cylindrical 
coordinate values x and 6 for the node, and vv is the node perturbation vector {vv(l), vv( 2 
vv(6)} that is used to update the location of the node. Variables vv(l), vv(2) and vv(3) correspond 
to perturbations to the axial, circumferential, and radial coordinates, respectively, and vv(4), 
vv(5) and vv(6) correspond to perturbations in the rotations about the axial, circumferential, and 
radial coordinate directions. During execution, the DIMP subroutine first reads in the measured 
imperfection data w 0 (x, 6) (e.g., data in Figure 3) from an external file and stores it in an array or 
look-up table. Then DIMP obtains finite-element nodal coordinate values x and 6 for a given 
node from the XYs array. A linear interpolation algorithm is used in DIMP to calculate the value 
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of the radial perturbation vv( 3 ) for each finite-element node based on the nodal coordinate values 
and the measured imperfection data stored in the look-up table. Similar calculations are also 
performed in DIMP to apply local rotations of the finite-element node that are associated with 
local changes in shell-wall curvature due to the imperfection, vv(4), vv(5). This process is 
repeated for each node called by the subroutine. The iunit, inode, prop , XYZg variables were not 
used. The actual code inside the subroutine is not reported here for brevity but can be found in 
Appendix D of Reference 26. 

Similarly, the shell-wall thickness t, mid-surface eccentricity ecz, and lamina material properties 
Ei, G 12 , V 12 and E 2 are adjusted at each integration point of each element in the finite-element 
models by using the user-written subroutine WALL. The WALL subroutine uses the following 
subroutine statement and common block for passing modeling information between the 
subroutine and the main STAGS program (variables used in the subroutine are highlighted in 
bold font): 


subroutine WALL( iunit, ielt, kelt, XYZg, XYs, zeta, ecz, ilin, iplas) 

COM M ON /W A LL 1/ lL, elL, ul2L, gL, e2L 

During execution, the WALL subroutine first reads in the measured thickness data t 0 (x, 6) (e.g., 
data in Ligure 4) from an external file and stores it in an array or look-up table. Then the WALL 
subroutine obtains finite-element integration point coordinate values x and 6 from the XYs array. 
A linear interpolation algorithm is used in WALL to calculate the individual ply thicknesses at 
each finite-element integration point based on the integration point coordinate values and the 
measured thickness data stored in the look-up table and are then adjusted using the tL vector (the 
size of the vector is based on the total number of plies, in this case, eight) The thickness of each 
lamina ply is assumed to be equal to the total measured laminate thickness divided by the total 
number of plies in the laminate, tL(i) = t 0 (x, 6)/S where i = 1, 8. The shell-wall mid-surface 
eccentricity ecz is calculated relative to the average measured shell-wall mid-surface as 
illustrated in Lig. 8; that is, ecz(x, 6) = -0.5(t ave - t 0 (x, d)). The updated lamina material 
properties Ei, G 12 , V 12 and E 2 are provided via the elL(i), u!2L(i), gL(i), and e2L(i) vectors, 
respectively, and calculated based on the rule of mixtures. In general, the rule-of-mixtures 
calculations used here assumed that any variation in the lamina ply thickness from the nominal 
thickness is due to a variation in resin volume only, and that the fiber volume remains constant 
for each ply. 

Several simplifying assumptions and approximations were made in modeling the geometry and 
stiffnesses of a lamina ply-gap detail. Lirst, the finite-element models are limited to modeling the 
shell-wall thickness variation as discrete changes, as illustrated in Lig. 8, and the resolution of 
the thickness variation is limited by the finite-element integration point spacing (e.g., 0.1 in.). 
Results from a study illustrating the potential effects of these modeling approximations and mesh 
refinement on the response of a shell with a ply-gap indicated that this modeling approach can 
affect the predicted response. As one might expect, these modeling approximations can lead to 
an artificial increase in the width of a lamina ply-gap, resulting in misrepresentation of the local 
bending stiffnesses and the mid- surface eccentricity of the shell wall at that particular point. In 
addition, the mesh refinement and integration point spacing used in the present models tend to 
provide models that are overly stiff in bending by a small amount when they are used for 
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modeling a ply-gap detail. However, these modeling approximations were shown to have a 
relatively small influence on the local bending response of the shell and only resulted in a 2% 
variation in the predicted buckling loads for the shells considered herein. 

The second assumption is related to the modeling of the stiffnesses of the lamina ply-gap detail. 
Two modeling approaches were considered. One approach is to model each ply-gap with a local 
reduction in thickness, as measured, including the local eccentricity ecz, and reducing the 
appropriate number of lamina plies in the shell-wall laminate model, hence, reducing the local 
stiffnesses associated with the local ply-gap detail. The second approach is to model the ply-gap 
with the local as-measured reduced shell-wall thickness and the corresponding local mid-surface 
eccentricity, but assume that the local thickness reduction is due to a reduction in the resin 
volume only and, consequently, keeping the fiber volume constant. This approach neglects the 
local stiffness reduction associated with the reduced number of plies. However, results from a 
study of this modeling assumption indicated that neglecting the local stiffness reduction 
associated with a ply-gap would cause only slight differences in the magnitude of the local 
bending response and no more than a 2% variation in the predicted buckling load for a 
compression-loaded shell. In addition, the results indicated that the local shell-wall mid-surface 
eccentricity is the most important feature of the ply-gap detail for these stability critical 
problems. As a result, the latter modeling approach was used in the present study. 

To provide a better simulation of the constraints provided by the potting material at the ends of 
the specimens, effective axial and radial potting- support stiffnesses were determined for each 
shell specimen using a two-dimensional generalized plane-strain finite-element analysis of the 
potting-material-shell-wall detail shown in Fig. 9a. Material properties of the potting compound 
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were characterized by Iosipescu tests of the potting material reported by Weiland et al. The 
nominal properties of the potting material are as follows: Young’s modulus E =1.15 Msi, shear 
modulus G 12 = 0.36 Msi, and Poisson’s ratio v = 0.59. The effective laminate axial stiffness E x 
was used in the shell-wall models which include E x equal to 11.83, 3.60 and 8.07 Msi, for shells 
Cl, C2 and C3, respectively. To determine the effective axial and radial potting-support 
stiffnesses, two numerical experiments were conducted using this local model. In the first 
numerical experiment, an axial displacement A was applied to the unpotted end of the shell 
model, and the resulting predicted axial strain response, shown in Fig. 9b, was used to calculate 
an effective axial stiffness Ka for the portions of the shell wall supported by the potting material. 
The predicted effective axial potted-shell stiffnesses are 26.55, 6.73 and 26.26 Msi for shells Cl, 
C2 and C3, respectively. In the second numerical experiment, a unit force F- was applied to the 
shell model, and the resulting w (radial) displacement and nodal force at the junction of the shell 
wall and the potting material was used to calculate an effective radial support stiffness Kr. The 
effective radial potting-support stiffness was predicted to be approximately equal to 1.0E5 lbf/in. 
The radial support stiffness was incorporated into the models via linear spring elements in 
STAGS. 

Nonuniform end loading of a specimen is attributed to initial specimen-end or loading-surface 
imperfections (see Figure 6) and to upper loading-platen rotations that are measured during the 
experiment. A typical set of measured loading platen rotations is shown in Fig. 7 and indicate 
that significant upper platen rotations occur from the onset of loading up to a load value of 
approximately 6,000 lbs. These rotations are attributed to an initial settling of the upper loading 
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platen onto the specimen. The rotations of the movable upper loading platen reach a steady state 
at a load value of 6,000 lbs., and the loading of the specimen, for the most part, continues 
without additional upper loading platen rotations from 6,000 lbs. up to the buckling load. During 
the buckling event, the upper loading platen undergoes an additional amount of rotation. 

These nonuniform displacements were simulated in the finite-element model in three parts. 
First, the measured upper and lower loading- surface imperfections d lop ( 6) and <5 bot(d), 
respectively, were included in the finite-element model by introducing an initial in-plane axial 
perturbation vvf 1) to each of the nodes at the loaded ends of the shell using the DIMP subroutine. 
Then, the compression load was applied to the shell in two parts. The nonuniform loading- 
surface imperfections, -d top ( 6) and -db ot ( 6), were applied as displacements to the upper and lower 
ends of the shell, respectively, in the first load step of the analysis to simulate a full contact 
condition between the shell ends and the loading platens. Then, the experimentally measured 
end-shortening displacement A and upper loading-platen rotations (p y and (p- were applied to the 
upper loading surface while holding the lower loading surface fixed as illustrated in Fig. 2; that 
is, u ( -L/2, 6) = A + R cos <py cost? + R cos <p z sint? -d top ( 6 ) and u (L/2, d) = -5b ot ( d) . 

These nonuniform loads were applied by using a user-defined subroutine USRLD in STAGS. 
The USRLD subroutine uses the following subroutine statement and call statement to define the 
applied displacement (variables used in the subroutine are highlighted in bold font): 

subroutine USRLD(iunit, X, Y, nrows, ncols, isys) 

call FORCE(P, It, Id, li, Ij, lax) 

where X and Y correspond to the x and 6 coordinates of the finite-element node and isys 
corresponds to the load set number. At the time of this work, the user was able to define two 
separate load sets in STAGS, isys = 1 was an active load set and isys = 2 was a passive load set. 
The subroutine FORCE is called inside USRLD to define the load for a particular node and 
includes the load magnitude P, load type It, and load direction Id. The load type can include 
applied displacement, applied point force, distributed load, live pressure, and thermal. 

Finally, the effects of selected specimen parametric uncertainties were determined based upon 
the results of a traditional combinatorial analysis of the effects of selected specimen parametric 
uncertainties. 17 The parametric uncertainties considered include uncertainties in the imperfection 
measurements, the lamina fiber volume fraction, and the applied load distribution. The 
imperfection measurement uncertainties are attributed to the accuracy tolerances of the 
coordinate measurement device used to measure the initial geometry of the shell. The shell-wall 
imperfection measurement values and the shell-end-shape imperfection measurement values 
used in the present study are accurate to within 0.0006 inches, or 0.75E-4R and 1.5E-4L, where R 
and L are the shell radius and length, respectively. The shell-wall thickness values are accurate to 
within 0.0012 inches, or 0.03 t nom ( t n0 m is the nominal laminate thickness equal to 0.04 inch). The 
lamina fiber volume fraction is specified by the manufacturer of the material to be equal to 0.65 
±0.03 for the 0.005-in.-thick graphite-epoxy preimpregnated tape material used to fabricate the 
test specimens. Applied load distribution uncertainties are measured indirectly by monitoring the 
measured and predicted axial strains at selected points near the top and bottom loading surfaces 
of the shell. A correction to the applied displacements can be determined from the differences in 
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the measured and predicted strains as follows. A user-written program external to the STAGS 
code was used to analyze the differences in the measured and predicted strains for a specified 
applied load value. This program used an iterative predictor-corrector method to determine a 
correction to the applied shell-end displacements. A new finite-element analysis was conducted 
with this displacement correction included in the model. This process was repeated iteratively 
until the difference in the measured and predicted strains reached a predetermined tolerance. A 
typical load correction is shown in Figure 10. 

Nonlinear Analysis Procedure 

The code uses both the modified and full Newton methods for its nonlinear solution algorithms, 
and accounts for large rotations in the shell by using a corotational algorithm at the element 
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level. The Riks pseudo arc-length path-following method - is used to continue a solution past the 
limit points of a nonlinear response. With this strategy, the incrementally applied loading 
parameter is replaced by an arc-length along the solution path, which is then used as the 
independent loading parameter. The arc-length increments are automatically adjusted by the 
program as a function of the solution behavior. The transient analysis option in STAGS uses 
proportional structural damping and an implicit numerical time-integration method developed by 
Park. - Additional information on the transient analysis procedure can be found in Reference 30 
by Riks, Rankin and Brogan. It should be noted that these advanced nonlinear solution 
algorithms (e.g., the corotational algorithm, Riks arc-length method, transient dynamic analysis, 
etc.) were first developed for STAGS and represented significant breakthroughs in the analysis 
for thin-walled shells at the time of their development, many of which have since been migrated 
to commercial finite-element codes. 

The prebuckling, buckling and postbuckling responses of the shells presented herein were 
determined using the following analysis procedure. A modified version of Riks pseudo arc- 
length path-following method in STAGS was used to compute the initial shell response until just 
before buckling occurred. The Riks subroutine was modified inside the STAGS sourcecode so 
that it could monitor the number of negative roots in the tangent stiffness matrix and make 
adjustments to the solution. More specifically, when a negative root was identified, indicating an 
unstable equilibrium solution, the solution was discarded and the analysis was automatically 
restarted from the previous stable equilibrium solution and would proceed with a reduced step 
size. This process would continue until a solution close the limit point was obtained based on a 
user-specified solution tolerance (typically less than 0.5% of the expected critical end- 
shortening). The unstable buckling response of the shell was predicted using the nonlinear 
transient analysis option of the STAGS code. The transient analysis was initiated from an 
equilibrium state close to the limit point by incrementing the end displacement by a small 
amount. An initial time step of 1.0E-8 seconds was used in the analysis, and the time step was 
automatically adjusted by the program as a function of the solution behavior. The transient 
analysis was continued until the kinetic energy in the shell had dissipated to a negligible level, 
which indicated that the transient response had attenuated. Once the transient analysis had 
attenuated to a near-steady- state solution, a quasi-static load relaxation option of the code was 
used to establish a stable postbuckling equilibrium state. 
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IV. Results and Discussion 


Numerically predicted and experimentally measured results for three compression-loaded 8-ply 
graphite-epoxy laminated cylindrical shells are presented in this section. The shell-wall laminates 
include two different orthotropic laminates and one quasi-isotropic laminate. The predicted 
results were obtained from finite-element models of geometrically perfect shells and shell 
models that include initial geometric shell-wall imperfections, shell-wall thickness variations and 
thickness-adjusted lamina properties, local shell-wall lamina ply-gaps, elastic boundary support 
conditions, and nonuniform loading effects, as described in section III, referred to herein as the 
imperfect shell. In addition, some models also include the effects of specimen parametric 
uncertainties. The results presented include predicted and measured load-end-shortening 
response curves, predicted and measured load-strain response curves, and predicted prebuckling, 
buckling and postbuckling deformation response patterns. 

Predicted and Measured Response Comparisons 

Analytically predicted and experimentally measured load-end- shortening response curves for 
shells Cl, C2 and C3 are presented in Fig. 11. Values of the axial load P and the end-shortening 
A are normali z ed by the linear bifurcation buckling load of quasi-isotropic shell C3, P^= 42,590 
lbs., and the nominal shell-wall thickness t nom = 0.04 in., respectively. The solid and dashed lines 
in the figure represent experimentally measured and analytically predicted results, respectively. 
Each shell has two predicted response curves representing analytically predicted upper and lower 
bounds to the imperfect shell response based on specimen parametric uncertainties described in 
section III. The measured results indicate that shells Cl, C2 and C3 exhibit global buckling at 
normalized load values of 0.652, 0.749 and 0.803, respectively, and are 7.8, 13.7 and 17.6% 
lower than the predicted loads for the corresponding geometrically perfect shells, respectively. 
The global buckling response of the shell is characterized by a sudden reduction in the axial load 
in the shell that is associated with the transient collapse. These results indicate that, for the most 
part, the measured load-end-shortening response curves for the shells fall on or within the 
predicted response bounds. 

The load-time history of the predicted transient collapse response of the imperfect shell C3 is 
shown in Fig. 12. For this shell model, buckling occurs at a normalized load of 0.890. The load- 
time history curve exhibits a sudden reduction in axial load until the collapse response attenuates 
and the axial load achieves a steady-state value. The kinetic energy in the shell obtains a 
maximum value during the transient collapse response and dissipates over time and the shell 
reaches a stable postbuckling equilibrium state after approximately 0.007-0.008 seconds. The 
shell exhibits postbuckling load-carrying capacity, however, the effective axial stiffness of the 
specimen is significantly reduced in the postbuckling range of loading. In addition, the results 
indicate that the postbuckling load level is relatively insensitive to the different types of 
imperfections considered herein. 

The corresponding transient deformation responses for selected time steps during the transient 
collapse response of the imperfect shell C3, indicated by the letters A through D in Fig. 12, are 
presented in Fig. 13a through 13d, respectively. Just before buckling occurs, the shell-wall 
deformations are characterized by several localized ellipse-like buckles located near one end of 
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the shell, as indicated in Fig. 13a. The localization in the deformation pattern is caused by the 
combination of a local geometric shell-wall imperfection and the bending boundary-layer 
deformations that form near the end of the shell. After approximately 0.0012 seconds have 
elapsed in the transient response, a single ellipse-like buckle has grown in amplitude and 
interacts with the destabilizing stresses in the shell wall to cause the general instability and 
collapse of the shell. The magnitude of the shell-wall radial displacement varies between ±0.5 t nom 
at this point in the response. After additional time has elapsed in the transient collapse response, 
additional local buckles have formed around the circumference and along the length of the shell, 
as indicated in Fig. 13c, and the normalized axial load has decreased from 0.890 to 0.554. In 
addition, some of the buckles in the shell begin to coalesce into larger diamond- shaped buckles. 
The magnitude of the shell-wall radial displacement varies between +3 t nom to -7 t nom . After 
approximately 0.01 seconds have elapsed in the transient response, the kinetic energy in the shell 
has dissipated to a negligible level indicating that the transient response has attenuated, and the 
shell has deformed into a stable postbuckling mode-shape, as indicated in Fig. 13d. 

Predicted initial post-collapse normal displacement contours and the corresponding observed 
moire fringe patterns for shell C2 (qualitatively similar to the observed response for C3) are 
shown in Fig. 14. These results indicate that the shell collapses into a diamond-shaped global 
buckling pattern with 14 half-waves around the circumference and two half-waves along the 
length, as predicted by the transient analysis. In addition, the predicted mode-shape is in-phase 
with the observed mode-shape. Similar results for shell C3 indicate that the shell collapses into a 
diamond-shaped pattern with 16 circumferential half-waves and two axial half-waves. However, 
the analytical results predicted that the deformation pattern can change its circumferential 
alignment by approximately 15° depending on the values of the parametric uncertainties used in 
the analysis. 

Analytically predicted and experimentally measured back-to-back strain response curves for a 
representative location on shell C3 are presented in Fig. 15. Values of the axial strain s x are 
normalized by the axial strain at buckling e cr = 0.00257. The strain results represent back-to- 
back surface strain measurements within the bending boundary layer near one end of the shell 
from strain gages 83 and 84. These gages are outer and inner surface strain gages located at x/L T 
= 0.071 and 6 = 180°. The bending boundary layer behavior is evident in the divergent character 
of the back-to-back predicted and measured strain results. These results indicate that the 
measured prebuckling strain responses, for the most part, fall on or within the predicted response 
bounds. However, the results also indicate that the predicted post-collapse strain response can be 
very sensitive to specimen parametric uncertainties. In addition, some discrepancy is expected 
based on sensitivity of the strain measurement in a region that exhibits large magnitude bending 
gradients. The measured and predicted prebuckling strain results for orthotropic shells Cl and C2 
show similar agreement (not shown here). 

V. High-fidelity Models and Analysis-Based Buckling Knockdown Factors 

Modern nonlinear analysis tools, advanced measurement and test technologies, and improved 
understanding of the effects of imperfections on the buckling of thin-walled shells have enabled 
the development of high-fidelity structural models that can predict the buckling response of these 
shells to a high degree of accuracy. Validated high-fidelity models, such as those presented in 
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section IV, are enabling the development of new analysis-based design knockdown factors for 
buckling-critical structures. Selected results that illustrate the use of high-fidelity models for the 
development of analysis-based knockdown factors are presented in this section. First, measured 
initial geometric imperfection data from six graphite-epoxy shells similar to shells C1-C3 are 
used to determine a manufacturing-process-specific imperfection signature for these shells. The 
imperfection signature represents a "first- approximation" imperfection shape that is suitable for 
developing preliminary-design data. This imperfection signature is then used as input into 
nonlinear finite-element analyses to predict the buckling response of a family of shells, based on 
shell C3, with different simulated imperfection shapes. These results presented in this section are 
documented in Reference 18. 

Manufacturing Imperfection Signature 

Shell-wall geometric imperfections that had been measured previously (see Refs. 15-17) for the 
six laminated-composite shells were used to derive a manufacturing-process-specific "first 
approximation" imperfection signature for this type of structure. First, a reference or "best-fit" 
cylinder for the initial imperfection measurement data for each shell was determined. 
Specifically, a data reduction program that uses the method of least squares was used to compute 
the rigid-body translations and rotations and the mean radius of each cylinder. Then, the 
measured radial imperfections, deviations from a perfect cylinder, were calculated with respect 
to the corresponding new best-fit cylinder. The measured imperfection data was then represented 
by a two-dimensional half-wave Fourier cosine series given by 

w 0 (x, 0) = t nom Y.k = o IfLo cos ^ cos sin InG ) (1) 

where R , L , and t nom are the shell radius, length and wall thickness; x and y are the axial and 
circumferential coordinates; Au and Bki are amplitudes of the different modal components in the 
Fourier series; and k and / are integers denoting the number of axial half-waves and 
circumferential fullwaves, respectively. Using a standard representation of this type for the 
measured imperfection enables convenient analysis and comparison of different imperfection 
distributions from different shells and different manufacturing processes. For example, a plot of 
the coefficient distribution of a half-wave Fourier cosine representation of the measured 
geometric imperfection shape for shell specimen C3 is presented in Fig. 16. The results indicate 
that the largest-magnitude component of the imperfection is the k = 0 and 1 = 2 component and 
that the magnitudes of the other components decrease rapidly as the wave numbers k and / 
increase. The k = 0 and 1 = 2 component corresponds to a cylinder with an elliptical cross- 
section that has a magnitude of 1.2 t nom . The next largest component is associated with k = 1 and / 
= 2 and has a magnitude of approximately 0.2 t nom . In general, the dominant coefficients are 
associated with the long wave-length imperfections in the circumferential and axial directions. 
Coefficients for values of k and / > 10 have a small influence on the overall imperfection shape, 
however, the short wave-length components are retained in the imperfection model since these 
modes may be associated with high imperfection sensitivity for the shells considered herein even 
though their magnitudes are generally very small. Results for shells Cl, C2, and C4-C6 indicate 
similar coefficient distributions and suggest that a manufacturing-process-specific imperfection 
signature for graphite-epoxy laminated-composite shells can be established in terms of a 
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characteristic Fourier coefficient distribution and range. To this end, a mean imperfection shape 
and a corresponding standard deviation to the imperfection shape were determined based on the 
measured imperfection data from the six shells, Cl - C6. The mathematical expectation or mean 
value of w 0 (x, 6) is given by 


p w (x,6) = t nom Y I k=oY.i=o cos — (.( A ki)coslne + (B kl )smlnd) ( 2 ) 

where (...) denotes the mean value. A contour plot of the normalized mean geometric 
imperfections ia w (x, 6 ) is shown in Fig. 17. The mean shell-wall imperfection u w is normalized 
by the nominal shell-wall thickness t nom = 0.04 inches. These results indicate that the initial 
geometric shell-wall imperfection is periodic in the circumferential direction, exhibits slight 
variations in the axial direction, and is similar in character to the measured imperfection shape 
for shell C3 shown in Fig. 3. The corresponding standard deviation of the geometric imperfection 
shape o w (x,6 ) is defined by 

o£(wo(*»0)) = <(w o (x,0) - < w o (x,0)» 2 > (3) 

A contour plot of the corresponding normalized standard deviation of the geometric imperfection 
shape is shown in Fig. 18. The results indicate that the maximum standard deviation is equal to 
2.172 t nom and is primarily associated with long wave-length components of the imperfection 
shape. 

Response of Compression-Loaded Shells 

The effects of four different imperfection shapes on the response of shell C3 are presented. 
Specifically, the effects of actual measured shell-wall imperfections, and three simulated 
imperfections that are based on the imperfection signature. The simulated imperfections include 
the mean imperfection and the mean imperfection plus or minus one standard deviation. A shell 
with the mean imperfection shape is denoted by p, w ; a shell with the mean imperfection shape 
plus one standard deviation is denoted by u w + o w ; and a shell with the mean imperfection shape 
minus one standard deviation is denoted by p w - o w . Load- shortening response curves and shell 
deformations are presented subsequently for the quasi-isotropic shell C3 that illustrate the overall 
response characteristics for the compression-loaded shell with the four different imperfection 
shapes. 

Four predicted load-shortening response curves are presented for shell C3 in Fig. 19 and include 
results for a shell with the actual measured imperfection shape;. The axial load P and end- 
shortening A are normalized with respect to the linear bifurcation buckling load of the 
geometrically perfect nominal shell, Puf = 42,590 lbf, and the nominal shell-wall thickness, t nom 
= 0.04 in, respectively. The load-end-shortening curves indicate a linear prebuckling response. 
The general instability occurs in the shell at normalized axial load levels that range from PI Pbif 
= 0.855, for a shell with u w + o w imperfection shape, to 0.943 for the shell with a u w - o w 
imperfection shape. 

Initial buckling deformations for shell C3 with the mean imperfection shape, the u w + o w 
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imperfection shape, and jx w - o w are presented in Figs. 20a through 20c and indicate similar 
response characteristics to those exhibited by the shell with the actual measured imperfection. In 
particular, the initial buckling deformation response of a shell with the mean imperfection and 
the u w + o w imperfection are similar to that exhibited by the shell with the actual measured 
imperfection and are characterized by a localization of a single buckle near one end of the shell. 
These results are significant because they suggest that the use of a mean imperfection shape and 
the corresponding standard deviation can bound the experimental buckling loads as well as 
predict physically meaningful response characteristics such as the buckling deformations. The 
correlation would improve still further if other measured imperfections were included in the 
finite-element model; e.g., measured thickness variation, material property variation, and loading 
nonuniformities, as shown in Refs. 17 and 18. Specifically, the measured thickness imperfection 
has been shown to account for an additional 5-10% reduction in the buckling load. In 
comparison, the lower-bound design recommendation given in NASA SP-8007 for the 
corresponding isotropic shell with R/t = 200 is equal to PJ P bi f = 0.47, which is very 
conservative, as expected. 


VI. Concluding Remarks 

The modem high-fidelity analysis approach offers a relatively affordable method to develop 
updated design knockdown factors for modern buckling-critical shell structures. The suggested 
design and high-fidelity analysis procedures presented herein should be used with a selected 
number of carefully conducted experiments that would be used to validate the design and 
analysis results. This approach can help avoid testing the large number of replicates of a design 
needed to develop empirical design factors. Furthermore, this nonlinear analysis procedure can 
be used as a parametric tool in the early stages of a design development program to determine 
the sensitivity of the response characteristics of a specific design to a number of different types 
of imperfections or differences in the idealized as-designed shell structure and the actual as- 
manufactured shell structure. Since this nonlinear analysis procedure can predict the local 
stresses and strains in a shell for any point in the load-response history of the shell, failure of the 
shell due to local stress or strain gradients could also be predicted. 

In recent years, a comprehensive project has been established at NASA to develop and validate 
new analysis-based shell buckling knockdown factors and design recommendations for modern 
launch vehicle structures. These new knockdown factors are based on the results from high- 
fidelity analysis tools such as the ones presented here and are being rigorously validated through 
sub-scale and full-scale structural tests. These new factors will account for the effects of 
geometric imperfections, loading nonuniformities, orthotropy, longitudinal welds, and combined 
mechanical, thermal and pressure loads typically found in large metallic cylindrical structures. It 
is expected that the new factors will enable significant mass reduction in the next generation of 
NASA launch vehicles as well as reduce development cost and risk. The STAGS finite-element 
code has played a prominent role in all phases of the development of high-fidelity analysis 
methods for buckling-critical structures and has enabled the development of these new analysis- 
based knockdown factors. It is with this in mind that one can recognize the continued need for 
research-oriented codes such as STAGS as a development test bed and benchmark for future 
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high-fidelity physics-based models that will be needed to meet the next generation of structural 
analysis challenges. 
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Figure 1. Typical laminated composite cylinder. 
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Figure 2. Finite-element model. 
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Figure 3. Measured geometric imperfection of quasi-isotropic shell C3. 
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Figure 4. Measured shell-wall thickness of quasi-isotropic shell C3. 
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Figure 5. Typical lamina ply-gap and outer-surface ridges in 8-ply laminated composite cylinders. 




Figure 6. Measured loading surface imperfection for shell C3. 
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Figure 7. Typical measured loading platen rotations during experimental testing. 
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Figure 8. Idealized ply-gap model detail. 
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Figure 9. Finite-element model of the potting-shell detail and typical predicted axial strain distribution at 
the mid-surface of the shell wall. 



Figure 10. Numerically predicted displacement corrections 8 c . top and 8 c _ bot for the top and bottom loading 
surface, respectively, inferred from experimentally measured axial strains for shell C3. 
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Figure 11. Predicted and measured load— end-shortening response curves for shells Cl, C2, and C3. 
Predicted results represent response bounds. 



Figure 12. Predicted load versus time history for the transient collapse of geometrically imperfect shell 
C3. 
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Figure 13. Predicted buckling response for an imperfect compression-loaded quasi-isotropic shell C3. 
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Figure 14. Observed and predicted initial post-buckling radial displacements for shell C2 (similar to the 
response for shell C3). 


25 

American Institute of Aeronautics and Astronautics 




Figure 15. Typical predicted and measured back-to-back strain response for quasi-isotropic shell C3. 
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Figure 16. Coefficient distribution of a Fourier series representation of the measured geometric 
imperfection for shell C3. 
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Figure 17. “First approximation” imperfection signature. 
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Figure 18. Standard deviation to the imperfection signature shown in Figure 17. 
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Figure 19. Predicted nonlinear response for the compression-loaded shell C3 with selected initial 
imperfection shapes. 



Figure 20. Predicted initial buckling deformation responses for the compression-loaded shell C3 with 
selected simulated initial imperfection shapes. 
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